Exactly Solvable Topological Chiral Spin Liquid with Random Exchange 
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We extend the Yao-Kivelson decorated honeycomb lattice Kitaev model [Phys. Rev. Lett. 99, 
247203 (2007)] of an exactly solvable chiral spin liquid by including disordered exchange couplings. 
We have determined the phase diagram of this system and found that disorder enlarges the region of 
the topological non-Abelian phase with finite Chern number. We study the energy level statistics as 
a function of disorder and other parameters in the Hamiltonian, and show that the phase transition 
between the non-Abelian and Abelian phases of the model at large disorder can be associated with 
pair annihilation of extended states at zero energy. Analogies to integer quantum Hall systems, 
topological Anderson insulators, and disordered topological Chern insulators are discussed. 
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I. INTRODUCTION 

Quantum spin liquids in two dimensions are often dif- 
ficult to firmly establish theoretically because of their 
strongly correlated nature^. An important advance in 
this direction was made by Kitaev who discovered a 
class of spin models with spin liquid ground states that 
can be solved exactlj*^. These models offer an existence 
proof for many proposed exotic spin liquid states, such 
as those with non-Abelian excitations^i^. A number of 
closely related models have been studied in two and three 
dimensions^, and controversial states such as such as two 
dimensional gaplcss spin liquids with a stable spin Fermi 
surface have been established^. Besides the existence of 
certain types of quantum spin liquids, their response to 
local impurities, and disorder more generally, is an im- 
portant and largely unexplored issue^.To date, only a few 
studies of local defects and impurities in Kitaev models 
have been publishecU^. 

In this article we examine an exactly solvable fully dis- 
ordered version of a Kitaev model, Eq.(IT]). on the dec- 
orated honeycomb lattice shown in FiglTJ As was orig- 
inally suggested by Kitae-\i^ such a lattice with triangu- 
lar plaquettes could lead to a chiral spin liquid ground 
state. In the clean limit, this model is known to pos- 
sess both Abelian and topological (finite Chern number) 
non-Abelian gapped chiral spin liquid phases^^— . In this 
work, wc focus on how random exchange disorder affects 
the phase boundaries and show there are analogs to the 
recently studied topological Anderson insulatorsii"— and 
disordered Chern insulatorsi^. Our main result is that 
disorder enlarges the parameter space of the topologi- 
cal phase and therefore can drive a transition into the 
topological phase. This implies that disorder could play 
a useful role in experimentally realizing certain classes 
of topological quantum spin liquids. We study in detail 
the Abelian and non-Abelian phases, and the transition 
between them as a function of disorder by examining (i) 
the spatial distribution of thermal currents under a small 
temperature gradient and (ii) the energy level statistics. 
We show that in the large disorder limit the transition 
between the two phases can be understood as the pair 



annihilation of extended bulk states at zero energy, anal- 
ogous to the integer quantum Hall effect (IQHE). To the 
best of our knowledge, a study of random disorder in a 
chiral Kitaev model has not been carried out before. 

The paper is organized as follows. In section II we 
introduce the Yao-Kivelson model, its mapping to non- 
interacting Majorana fermions and its different ground 
states. Next, we introduce the form of exchange dis- 
order that we have studied. In section III wc discuss 
how the Chern number of the disordered system is de- 
termined through the numerical calculation of the linear 
response thermal Hall conductivity. We then present the 
disordered phase diagram of this system. Next we vi- 
sualized the local thermal currents in the topologically 
non-trivial phase at large disorder. In section IV we re- 
port the results of a energy level statistics analysis that is 
used to characterize the disordered Majorana mode wave- 
functions at various points in the phase diagram. Finally, 
in section V we comment on the behavior of the spectral 
gap and thermal conductivity at large disorder and near 
the quantum phase transition. 



II. MODEL 

We consider spin- 1/2 moments on the lattice shown in 
FiglT] with Hamiltonian, 
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where /i = x,y,z,x' ,y' ,z' denotes the various types of 
links in the lattice, labelled as in FigH] The J,J'>0 are 
exchange couplings between the spins represented by the 
Pauli matricies a"', a = x,y,z on site j. The presence 
of plaquettes with an odd number of sides results in the 
spontaneous breaking of time-reversal symmetry^i^. 

The Hamiltonian ([T]) is solved by introducing Majo- 
rana operators Cf j 'Cf? 4" j Ci at each site »^. The combina- 
tion a^ = i^^Ci faithfully reproduces the Pauli operator 




FIG. 1. (a) The decorated honeycomb lattice and the labehng 
of hnks describing the interactions in Eq.(IT]). (b) The Z2 
gauge field Uij configuration and their ground state fluxes^. 
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where A = S,f^^£,fci. Recasting ^ 



in terms of Majoranas, yields a "tight-binding" Hamilto- 
nian for c in a static background Z2 gauge field 

{ij)GAi-links (ij)GM'-links 

(2) 
where Uij = — i^f^^ with fi ~ x,y, z and fi' = x',y',z'. 

PHP = H then recovers ([T]). The Uy- commute with H 
and thus provide good quantum numbers Uij = —Uji = 
±1. They are non-dynamical Z2 gauge fields. Thus, their 
fluxes are gauge invariant and integral invariants of mo- 
tion of H. The H ground state (which is degenerate) is 
specified by the flux configuration which has the least en- 
ergy when no excitations are present. Yao and Kivelsor*^ 
have shown that the minimal flux configuration is the 
one illustrated in FiglT] (b). In the clean system, the 
ground state of H for J' /J < -^3 is gapped and topolog- 
ically non-trivial with Chern number t/ = ±1, while for 
J'/ J > \/3 the ground state is gapped with v — 0^. 

It was discussed in refi^ that within a flux sector, P 
projects out half the possible eigenstates of H and this is 
believed to be a general feature of all Kitaev models^— li^. 
In particular, the physical states will depend on the rel- 
ative sign of the exchange couplings J and J', the global 
boundary conditions {e.g. on a torus), the specific flux 
sector and whether the number of c-field excitations is 
even or odd. Wc have verified our ground states are 
physical by numerically applying the projection proto- 
col of Refi^ to the Yao-Kivelson model for finite system 
sizes. We then considered a disordered version of ([T]) by 
allowing J and J' to vary randomly from link to link as 
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<ij> 



J (1 + X<ij>) , J^y>, = / (1 + X<ij>0 , (3) 



where < ij > are any of the x, y, z-links and < ij >' any 



of the 



,y',2' 



■links. The random variables X 



<ij> 



have 



uniform distribution P (|X<ij>| < W/2) = 1/W and zero 
otherwise. The parameters j and j' set the average; W is 
dimensionless and is limited to M^ < 2 so that J, J' > 0. 
By adiabatic continuity, this guarantees ^ remains in 
the ground state flux sector shown in Fig[T] (b). 



III. THERMAL CONDUCTIVITY 

We consider a finite rectangular system with 18 x 27 
unit cells (each with with 6 spins) with N = 2196 sites 
and open boundary conditions. For a finite system with 
h' = ±1 for some choice oi j,j' and disorder realization, 
an edge mode should be present at zero energy in the 
bulk gap>^. In contrast to the IQHE, these states do not 
couple to conventional electromagnetic fields, but never- 
theless do coupled to thermal gradients. The topological 
phase with z/ = ±1 would be manifested as a quantized 
thermal Hall conductivity, a^^y, with a^y = ^^kl^Tv, 
where fc^ is Boltzmann's constant and T the absolute 
temperature^. We computed cr*^ numerically for specific 
disorder realizations and then averaged over 100 realiza- 
tions. 
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FIG. 2. (Color online) A contour plot of (a^y)dis in units of 
^i-k%T. Finite size effects lead to a rounding and a small shift 
in the critical j'^^^^ = vSj expected in the clean limit for the 
transition between topological and non-topological phases^. 

Since Q with ([3]) is exactly diagonalizable, the con- 
ductivity can be determined within the Kcldysh Green's 
function formalism which has been successfully applied 



to mcsoscopic transport in quantum wires 
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We 



implemented this technique by coupling our disordered 
system to two ideal "leads" which are modeled as clean 
versions of the system with j' = j such that chiral edge 
modes are supported in the leads. If the disordered sys- 
tem under study possesses chiral edge modes, these and 
only these states will couple perfectly (without reflection) 
with the modes already present in the leads at low ener- 
gies. Thus a calculation to linear response of the thermal 
conductance, will then yield the quantized conductance. 
If however the system does not possess chiral edge modes 
and has i/ = which is different form that of the ideal 
leads, the thermal transmission will be heavily or entirely 
suppressed. Wc should also point out that if the leads do 
not support low energy chiral edge modes, then the thcr- 
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FIG. 3. (Color online) Local thermal current density profiles, /*]\ of a disordered quasi one-dimensional "wire" (running right 
to left) for typical realizations possessing chiral edge modes at large disorder W = 1.75. (a) j' = 1.5j (b) j' — 2.1j. The currents 
are displayed in units of ^j-k%T. Currents below 0.5 are shown as translucent and vary from fully transparent (/* = 0) to 
fully opaque (llf = 0.5). Only a part of the system is shown in fine detail, but the respective insets contain of coarse grained 
current profiles of the entire "wire" lattice between the two "leads" . 



mal conductance will always be trivially zero. One could 
also use a different system altogether to model the leads, 
say the honeycomb Kitaev model in its gapped chiral 
phase. The crucial requirement is that it is gapped in 
the bulk but has edge modes of the same number and 
chirality as the modes of the clean system in the non- 
Abelian phase. 

The right lead is then set to T = and the left is kept 
at a very small T > 0. To linear order in T, cr*^ is given 
by a Chester-Thellung-Kubo-Greenwood ioiimila^^^^—, 
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where the single-particle retarded and advanced Green's 
functions are evaluated at zero energy (E = 0) and the 
trace is taken over transverse modes. The leads, which 
serve as a heat source and sink, are modeled by self 
energy matrices, 5]('''^*/'''s'i')'(''°*/'^'^^\ which only act on 
sites that interface with the leads. These extra terms 
are included in the calculation of G™*/'''^^, and T^<^f^/"si^^ 
are associated with the imaginary parts of 5]('°^'/'''sht),rct^ 
The results of these calculations yielded the conductivity 
contour plot in Fig[2l which also serves as a phase dia- 
gram. The most striking feature is the expansion of the 
topologically non-trivial phase region for intermediate to 
strong values of disorder strength, 1 < W < 2. Similar 



effects have also been observed in disordered topological 
insulatorsiir— and disordered Chern insulatorsii. 

We visualized the local thermal current density across 
the Hnk < ij >, J'*^, in Fig. [3] for large W disorder 



±1 



realizations in the non-trivial topological phase v 
by adapting the methods of Reff^^. The If^ are deter- 
mined from the non-equilibrium Green's function G^ in 
the linear response regime of small T. Figini^a) clearly 
demonstrates the chiral nature of the edge modes which 
are localized on the upper edge. However, as j' is tuned 
towards criticality as in Fig|3l^b), the edge currents begin 
to impinge into the bulk. Eventually a backscattering 
channel opens when the currents connect the different 
edges. In close analogy to the plateau transitions of the 
IQHE, the topological phase transition can be associated 
with the closure of the bulk gap and the appearance of 
an extended bulk eigenstate at £■ = that percolates the 
system. 



IV. ENERGY LEVEL STATISTICS 

The results of the previous section suggest that even 
at large finite disorder, the phase transition in the bulk 
proceeds from a thermal insulator to a thermal insula- 
tor via a bulk thermal metal at criticality. This critical 
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FIG. 4. (Color online) Energy level statistics at various points in the phase diagram, (a)-(l) The variance of the energy spacing 
s shown as the solid blue (dark gray) line in units of mean level spacing 5 against the GUE expectation of 0.178 shown as the 
horizontal dashed green line. Also shown is the averaged density of states which is the solid red (light gray) line. Statistics 
are derived from 10'^ disorder realizations and binned in energy windows of width A_B — 0.25 j. Shown in the center are the 
locations of (a)-(l) in an approximate phase diagram. Spacing distribution histograms are taken from energy windows (indicated 
by the black lines and arrows) of special interest with comparisons to the Wigner GUE and Poisson distributions. The finite 
gap at _B = does not close in (a)-(c) which is not apparent here with this scale and bin size. 



metallic phase is manifested by the emergence of an ex- 
tended Majorana mode in the bulk at -B = 0. Thus with 
the aim of determining the localized/extended character 
of the Majorana mode wavefunctions, we have performed 
an energy level statistics analysis with comparisons to the 
expectations of random matrix theory. More specifically, 
a statistical analysis of closest energy level spacings, s, 
allows for an indirect determination of the spatial char- 
acter of the wavefunctions across the spectru m^"^'^^ . In 
addition an analysis based on level spectra has the tech- 
nical advantage of not requiring the explicit computation 
of disordered wavefunctions. Hence it has allowed us to 
efficiently study the entire spectrum for many points on 
the phase diagram without the need for intensive com- 
puting resources. Since H possesses particle-hole but no 
other symmetries, it falls into the ZJ-symmctry class of 
the ten-fold Altland and Zirnbaucr classification^^, which 
is also shared by Bogoliubov-de Gennes (BdG) Hamilto- 
nians with the same symmetries. Moreover, it is pre- 
dicted that at energies E much greater than the mean 
level spacing, 5, the energy level statistics will reduce to 
that of the Gaussian Unitary ensemble (GUE) of Wigner 
and DysonSS,. Delocalizcd eigenstates of similar energies 



are then expected to experience level repulsion described 
by the GUE Wigner surmise, P{s) = ■||s^e~^*"/'^, where 
s is measured in units of 6. Level repulsion is exempli- 
fied by the decay of probability density as s ^^ 0+. In 
contrast, strongly localized eigenstates of similar energies 
experience less level repulsion and are more likely to be 
described by Poisson statistics P{s) = e^^ By gathering 
statistics of level spacings from energy windows across 
the spectrum and making comparisons to these limits, 
one can distinguish between regions of strong localiza- 
tion, intermediate localization, and delocalization. We 
have performed such an analysis, and a summary of our 
results is presented in FiglH Spectra from disordered 
samples for systems of the dimensions mentioned earlier 
with periodic boundary conditions at various values of 
j' and disorder strength W were averaged over 10'^ re- 
alizations. Shown in Fig|4l[a-1) are the variances of the 
normalized energy level spacings Var(s) and the averaged 
density of states (DOS). 

For small disorder W = 0.25 the bands are fairly well 
defined by clear energy gaps. Within these bands lie re- 
gions of extended states distinguished by plateaus where 
the spacing variances agrees well with the GUE value 



0.178. Samples taken from these low regions of variance 
also indicate a strong agreement with the GUE distri- 
bution. Towards the band edges however, the distribu- 
tion of spacings tend to resemble Poissonian statistics, 
and more crucially the level repulsion is greatly reduced. 
With increasing W, the energy bands broaden and close 
some gaps in the DOS but there still remains gencrically 
a strong repulsion from E = except at a phase tran- 
sition. As W increases, the spacing variances across the 
spectrum increase and exhibit less level repulsion. Inter- 
estingly, the DOS and spacing variances on either side 
of the critical line behave rather differently. Specifically, 
the topological v = ±1 phase develops a narrow robust 
region of extended-like states with strong level repulsion 
located near E = which is marked by a clear dip in 
the spacings. As shown by the histogram taken from 
a bin in Fig|4l[b), states in this energy window exhibit 
strong level repulsion which matches the GUE distribu- 
tion fairly well. By increasing j' in Figs|4lja-d), we also 
observed that these extended states eventually merge or 
annihilate at the phase transition. These results then 
allow us to connect with the findings of the previous sec- 
tion. Namely, at the phase transition there emerges an 
extended E = state (s) which we associate with the 
pair annihilation at £' = of two extended regions which 
evolve from the bulk bands. We note that an identical 
phenomena of pair annihilation had also been previously 
observed in topological insulator, IQHE and topological 
superconductor systemsi^i^i. This system, however, is 
the first example of a disordered Kitaev model shown to 
exhibit this sort of physics. We remark that FigUl^d) in- 
dicates an intermediate distribution near the transitions, 
and Fig|4ljk) shows an interesting bimodal distribution 
in the vicinity of strong DOS fiuctuations. 



V. LARGE DISORDER AND FINITE SIZE 
EFFECTS 

In this last section we briefly comment upon the spec- 
tral and transport properties of the most interesting part 
of the phase diagram. That is at large disorder and across 
the phase boundary. We also discuss related finite size 
effects that are observed. 

First we study the gradual closing of the gap at large 
constant disorder strength. An enlarged plot of the disor- 
der averaged density of states (DOS) near E = around 
the transition is shown in FiglDJa). Spectra was taken 
from a system of size N = 2916 and disorder strength 
W = 1.75 with j' varying between 2.Qj and 3.0j. For 
these parameters and to within our numerical precision, 
the gap is seen to close significantly in a window between 
j' = 2.25J and j' = 2.5Qj. For comparison, the disorder 
averaged thermal Hall conductivities were computed for 
the relevant range of parameters and plotted in Fig 5(b). 
The conductivities of larger system sizes was also studied 
and is also shown in FiglSKb). The additional data allows 
us to put the critical j' at approximately 2.31j. This is 
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FIG. 5. (Color online) Spectral and transport properties in 
the vicinity of the phase transition at large disorder strength 
W = 1.75. (a) Enlarged plots of the disorder averaged den- 
sity of states (DOS) near i? = where j' varies between 
the topologically non-trivial phase (2.0j) to the topologically 
trivial (3.0j) for a system of fixed size A'^ = 2916. Each point 
was derived from 10'^ realizations, (b) The disorder averaged 
thermal Hall conductances < cr^^ > across the phase transi- 
tion and for various system sizes A'" but constant aspect ratio. 
Conductances are expressed in the quantized thermal conduc- 
tance unit i^kgT and were averaged over 500 realizations. 
The critical j' is inferred to be j'^rit ~ 2.31j 



the point where all the conductivity curves cross. 

Next we examined the DOS near i? = for various 
system sizes near and away from the transition. Shown 
in Figini are DOS plots near E = Q for system sizes 
N = 1296, 2916 and 5184. We remind the reader that 
the system size used in our main results is N = 2916 
sites and we have not spectrally studied larger sizes in as 
much detail due to current hardware limitations. Inter- 
estingly the system size dependence is affected by how 
close j' is to its critical value. Deep in the topological 
1/ = ±1 phase with j' = 1.5j as shown in FiglHl^a), the 
DOS varies weakly with system size near the gap. How- 
ever in Figini^b), when j' = 2.25 j and where the gap 
is seen to close, the average density of states at ii^ = 
diminishes noticeably with increasing system size. It sug- 
gests that the gap's closure at j' = 2.25j is really a finite 
size effect. Although, we believe that a more detailed 
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FIG. 6. (Color online) Disorder averaged density of states 
(DOS) near E — for various system sizes and large disorder 
strength W = 1.75. (a) Deep in the topological u — ±1 phase 
j' — 1.5j where a robust gap develops which is size insensitive, 
(b) In proximity of the phase transition j' = 2.25J where the 
gap closes but DOS(i? = 0) diminishes with increasing system 



spectral analysis with larger samples and sizes is required 
to firmly establish this, which is planned for future work. 
Thus we have confirmed that for a finite system the 
gap is closed at the transition. However, to within our 
numerical precision it remains closed over a finite window 



of parameters. But this window lies inside a region of 
rapidly changing conductivity which contains the critical 
j'. Furthermore, we believe this behavior of the gap to 
be a finite size effect. 

Finally we remark that while it is certainly necessary 
that the gap closes at the transition, it may not be 
a sufficient condition. In a disordered system where 
disorder broadening may occur, the gap may be closed 
by states that are localized but are not topological. In an 
IQHE system, this would correspond to the Fermi level 
being placed in a band of localized states. Nevertheless, 
our numerics indicate that at least for our choice of 
model disorder and at large disorder strengths, the 
E=0 gap remains robust in the topological phase and 
gradually closes at the transition, but not sharply. 

VI. CONCLUSION 



We report a disordered generalization of the exactly 
solvable Yao-Kivelson-Kitaev model which supports chi- 
ral spin-liquid ground states with both trivial z/ = and 
non-trivial i/ = ±1 topological phases. We mapped out 
a disorder phase diagram of this system and showed that 
exchange disorder can enlarge the region of the topolog- 
ical phase. We visualized the local thermal current den- 
sity of the topological phase in contact with heat reser- 
voirs. A statistical analysis of the energy level spacings 
at various regions of the phase diagram yielded evidence 
for the nature of the wavefunctions in agreement with the 
thermal current studies. The low energy gap was shown 
to be robust in the non-Abelian phase even at large dis- 
order and closes at the onset of the phase transition. Our 
work shows there are many interesting parallels between 
topological quantum spin liquids and other topological 
phases2^, and we expect more to be discovered in the 
future. 
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